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LONG-TERM  GOALS 

The  long-term  goal  of  this  research  is  to  facilitate  physically  motivated  parameterization  of  small-scale 
physics  in  practical  ocean  forecast  models.  It  is  envisioned  that  broad  application  of  robust,  easy-to- 
use  simulation  tools  in  which  non-hydrostatic  physics  is  directly  resolved  will  provide  a  framework 
within  which  parameterization  schemes  can  be  formulated  and  quantitatively  evaluated. 

OBJECTIVES 

The  primary  aim  of  this  proposal  is  to  refine  an  existing  set  of  numerical  tools  for  the  simulation  of 
non-hydrostatic  flow  over  rough  topography,  to  undertake  a  series  of  validation  and  demonstration 
simulations  using  these  tools,  and  to  make  the  tools  themselves,  their  configuration  scripts  and 
documentation  available  to  the  community  at  large  via  the  web. 

APPROACH 

This  modeling  effort  employs  efficient,  high-order  numerical  techniques  on  parallel  computers  to  solve 
the  non-hydrostatic  equations  of  motion  for  stratified  flow  in  the  presence  of  irregular  topography.  It 
represents  a  merger  and  compromise  between  two  separate  lines  of  code  development.  The  first  is 
based  on  spectral  methods  for  turbulent  stratified  flows  in  the  absence  of  boundary  effects  while  the 
second  uses  curvilinear,  boundary  fitting  coordinates  along  with  implicit,  compact  differentiation 
schemes.  In  this  effort,  the  underlying  numerical  grid  is  Cartesian  with  the  boundaries  represented 
independently  using  immersed  boundary  techniques.  A  configurable  mix  between  spectral  and 
compact  differentiation  is  used  in  combination  with  fourth  order  accurate  time  stepping. 

Harvey  Seim  (University  of  North  Carolina),  Don  Slinn  (University  of  Florida),  Pascale  Lelong  and 
Bob  Robins  (Northwest  Research  Associates,  Bellevue  WA),  Bill  Smyth  (Oregon  State  University), 
and  Miles  Sundermeyer  (University  of  Massachussetts,  Dartmouth)  have  all  contributed  substantially 
to  this  effort  by  applying  different  versions  of  the  models  to  their  application  problems  and 
communicating  strengths,  weaknesses  and  problems  back  to  me. 

WORK  COMPLETED 

Several  tasks  of  a  technical  nature  have  been  accomplished  in  the  first  fiscal  year  of  this  project.  A 
brief  summary  is  given  below. 
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1.  Generalized  representation  of  topography:  Early  implementations  of  the  code  assumed  that  the 
topographic  height  varied  only  in  one  spatial  coordinate.  This  was  a  necessary  simplification  for  the 
boundary  fitting  grid  approach  where  a  curvilinear,  orthogonal  mapping  between  spatial  and 
computational  coordinates  was  required.  It  was  introduced  purely  for  programming  convenience  and 
for  ease  of  testing  in  the  immersed  boundary  code.  Specifying  the  bottom  topography  is  now 
accomplished  within  a  single  user-defined  subroutine  that  prescribes  the  height  of  the  topography 
h(x,y),  relative  to  some  reference  level,  given  an  arbitrary  horizontal  location  (x,y).  The  sampling 
requirements  for  the  topographic  representation  are  internally  calculated;  the  user  need  only  specify  a 
continuous  form.  Similarly,  the  domain  decomposition  and  the  consequent  mapping  of  the  immersed 
boundary  points  to  different  processors  is  internally  generated  and  shielded  from  the  user. 

2.  Spectral-differentiation:  User  configurable  switches  are  now  in  place  to  permit  spectral 
differentiation  whenever  appropriate.  The  spectral  routines  from  the  original  codes  have  been  upgraded 
to  use  the  FFTW  routines  from  MIT,  which  have  been  tuned  for  a  wide  variety  of  architectures  and  are 
generally  the  fastest  portable  routines  available.  Using  immersed  boundary  techniques,  it  is  often 
possible  to  define  the  computational  domain  in  such  a  way  that  “dead- zones”  or  regions  of  now  flow 
can  be  tied  together  with  periodicity  conditions  while  the  interior  flow  of  interest  exhibits  no  such 
structure.  This  permits  the  accuracy  of  spectral  differentiation  schemes  to  be  employed  for  seemingly 
non-periodic  flows. 

3.  Direct  solver  for  pressure:  Generally,  a  three-dimensional  elliptic  equation  needs  to  be  solved 
for  pressure  at  each  time  step  in  a  non-hydrostatic  code.  This  portion  of  the  algorithm  is  both  expensive 
and  demanding  in  the  sense  that  errors  and  inaccuracies  can  lead  to  unacceptable  degradation  of  the 
overall  solution.  I  have  coded  up  a  parallel  implementation  of  a  multi-grid  iterative  solver  for  the 
generalized  elliptic  equation  required.  To  date,  this  approach  has  proven  unsatisfactory  for  problems 
characterized  by  large  spatial  aspect  ratios.  For  these  problems,  convergence  occurs  either  too  slowly 
to  be  practical  or  does  not  occur  at  all.  For  problems  amenable  to  spectral  methods  in  the  horizontal 
however,  a  far  more  efficient  solution  is  available.  For  these  problems,  the  horizontal  portion  of  the 
problem  can  be  solved  semi-analytically,  i.e.  using  discrete  Fourier  transformations,  leaving  only  a 
one-dimensional  problem  that  needs  to  be  solved  numerically.  It  is  feasible  and  efficient  to  solve  this 
Id  problem  directly,  i.e.  by  forming  and  factoring  the  matrix  explicitly  outside  of  the  main  time 
stepping  loop.  This  approach  has  now  been  coded  and  tested  and  is  implemented  automatically  if  the 
appropriate  periodic  structure  of  the  computational  domain  is  detected. 

4.  Port  to  several  machines:  The  codes  have  been  ported,  tested  and  run  on  several  different 
parallel  architectures  including  three  different  Cray  T3E  environments,  an  Athlon  Beowulf  cluster  with 
Myrinet,  a  PHI  Beowulf  cluster  with  Myrinet,  an  Alpha  cluster  with  Ethernet  and  serial  workstations 
running  Finux. 

5.  Adiabatic  immersed  boundary  conditions:  There  are  many  approaches  to  handling  the  boundary 
conditions  prescribed  on  the  immersed  boundary.  Generally,  a  subsidiary  problem  is  solved  for  a  field 
of  forces  that  are  to  be  applied  to  the  fluid  at  grid  points  in  the  vicinity  of  the  boundary.  The  strength  of 
these  forces  is  calculated  so  that  when  applied,  the  fluid  is  brought  to  rest  at  the  immersed  boundary.  A 
similar  approach  has  been  developed  for  handling  the  adiabatic  condition  for  scalars.  A  set  of  control 
points,  just  below  the  immersed  boundary  is  defined.  These  points  represent  locations  for  sources  or 
sinks  of  heat  with  prescribed  spatial  decay  characteristics,  i.e.  a  localized  hot  spot  that  decays  like  a 
Gaussian  in  each  space  dimension.  At  each  time  step,  the  values  of  the  scalar  field  at  the  grid  points  are 
interpolated  to  the  immersed  boundary  points.  These  values  are  then  used  to  determine  approximately. 
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whether  the  temperature  values  at  the  eontrol  points  should  be  inereased  or  deereased  to  guide  the 
solution  toward  satisfaction  of  the  boundary  conditions. 

6.  Parallelized  LES  closure:  The  present  implementation  is  based  on  strait- forward  Laplacian 
dissipation  and  diffusion.  An  earlier  version  had  the  option  of  a  stratified-form  of  the  Smagarinsky 
LES  closure.  Work  is  ongoing  but  not  completed  to  provide  a  wider  variety  of  SGS  closure  options. 

RESULTS 

A  version  of  the  immersed  boundary  method  was  combined  with  spectral  and  6'^  order  compact 
differentiation  schemes  in  an  algorithm  for  stratified,  rotating  flow  in  the  presence  of  irregular 
application  of  the  boundary  conditions  at  the  topography  z=h(x,y)  as  opposed  to  applying  a  linearized 
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Figure  1:  Instantaneous  view  of  generation  of  internal  tide  at  horizontal  scale 
27i/ko=2500m  at  critical  slope.  The  bottom  left  panel  shows  the  phase  of  the 
forcing  at  the  time  of  the  snapshot. 


condition  at  z=0  appears  to  be  important  for  finite  amplitude  topography. 

Eigure  (1)  shows  a  snapshot  of  various  correlations  in  an  x-z  plane  for  oscillatory  flow  over  a 
corrugated  bottom  at  critical  slope  (wave  angle  matches  topographic  slope).  Several  differences  are 
apparent  in  comparison  with  recent  analytical  work  on  internal  tide  generation.  In  particular, 
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Figure  (2)  illustrates  the  implementation  of  the  immersed  boundary  approach  for  an  adiabatic  condition 
on  density.  Figure  (2a)  shows  the  density  perturbations  in  an  x-z  plane  along  with  an  immersed 
boundary  defining  a  channel  extending  into  the  plane  of  the  figure.  The  snapshot  is  taken  shortly  after 
startup  in  a  simulation  where  an  along-channel,  geostrophically  balanced  flow  adjusts  to  the  presence 
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Figure  2.  Upper  panel  (a)  shows  the  density  perturbation  from  a  linear  profile, 
emphasizing  the  heat  source  just  below  the  immersed  boundary.  Lower  panel  (b) 
shows  the  total  density  field  that  satisfies  the  required  adiabatic  condition.  The 

vertical  coordinate  is  in  meters. 


of  no-slip  channel  walls  by  establishing  a  cross-channel  circulation.  At  the  time  of  the  snapshot,  the 
cross-channel  flow  is  weak  and  shown  only  schematically  by  the  arrow.  The  initial  density  field  was 
consisted  of  a  linear  profile  in  z,  thus  violating  the  required  condition  of  no  normal  gradient  at  a 
curving  boundary.  The  immersed  boundary  treatment  detects  this  flaw  and  introduces  heat  sources 
below  the  immersed  boundary  such  that  the  density  boundary  conditions  are  quickly  and  continuously 
maintained.  Figure  2(b)  shows  the  total  density  field  satisfying  the  adiabatic  boundary  condition.  (Note 
the  large  horizontal  to  vertical  aspect  ratio.  Though  the  right-to-left  near  bottom  flow  is  only  just  being 
spun  up,  its  presence  is  detectable  as  an  asymmetry  in  the  density  contours  just  above  the  boundary.) 

The  region  below  (or  to  the  outside)  the  immersed  boundary  constitutes  a  “dead- zone”  of  no  flow.  This 
permits  the  use  of  periodic  computational  boundary  conditions  and  hence  spectral  differentiation  and  a 
direct  solution  method  for  pressure. 
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IMPACT/APPLICATIONS 


It  is  the  intent  of  this  work  that  high-resolution,  non-hydrostatic  process  modeling  will  become 
accessible  to  those  analyzing  data  and  or  running  larger-scale  ocean  models  and  will  aid  the  ongoing 
development  of  physically-based  parameterization  schemes  for  these  models. 

TRANSITIONS 

Several  groups  are  using  the  codes,  at  different  stages  of  development  for  their  own  scientific 
applications.  In  addition,  the  NASA  Office  of  Earth  Science  has  funded  a  three-year  grant  entitled 
Ocean  Mixing  from  Shear  Deformation  of  Ice  Leads:  An  Important  Component  of  the  Surface  Heat 
Budget  of  the  Polar  Oceans  (PI  Max  Coon,  NWRA)  which  will  use  this  model  on  the  T3E  at  ARSC. 
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